//#include <dynamic_matrix.hpp>
#include <matrix.hpp>

#include <iostream>
#include <iomanip>
#include <algorithm>
#include <iterator>

#include <boost/lambda/lambda.hpp>

#include <random.hpp>

using namespace boost::lambda;
using namespace std;
using namespace sm;

int main()
{
    cout << setprecision(5);

    dynamic_matrix<long double> a(11,11);

    cout << a << endl;

    long double b[]= { 1, 2, 3, 4, 5, 6, 7, 8, 9 };

    copy( b, b+3, a.col_rbegin(1));

    copy( b+3, b+8, a.row_rbegin(4));
    
    dynamic_matrix<long double> c = a;
    fill( c.rbegin(), c.rend(), 100.0 );
    copy( b+3, b+8, c.col_rbegin(4));

    dynamic_matrix<long double> d = -c;
    

    copy( d.upper_diag_rbegin(4), d.upper_diag_rend(4), ostream_iterator<long double>(cout, "\t"));
    cout << endl;
    copy( d.lower_diag_rbegin(0), d.lower_diag_rend(0), ostream_iterator<long double>(cout, "\t"));
    cout << endl;
   
    c[5][5]= 13.3;

    cout << c << endl;
    cout << endl;

    dynamic_matrix<int> e = c;
    cout << e << endl;
    cout << endl;
    
    d+=c;
    cout << d << endl;
    cout << endl;
   
    dynamic_matrix<long double> f = d + 10.0;

    cout << f << endl;

    dynamic_matrix<long double> g = f*d;

    cout << g << endl;
    cout << endl;

    dynamic_matrix<long double> h = g/87.0;

    h += 1;

    cout << h << endl;
    cout << endl;
    
    h.transpose();

    h -=1;
    
    cout << h << endl;


    dynamic_matrix<long double> i = h/g;

    fill( i.diag_rbegin(), i.diag_rend(), 1.3 );

    dynamic_matrix<long double> j = 1.0 / i;

    cout << i << endl;
    cout << endl;
    
    cout << j << endl;
    cout << endl;

    dynamic_matrix<long double> k(17,17);
    fill( k.rbegin(), k.rend(), 0 );

    dynamic_matrix<long double>l(k);
    fill( k.diag_rbegin(), k.diag_rend(), 10 );
    fill( l.diag_rbegin(), l.diag_rend(), 0.1 );

    cout << l/k << endl;
  
    k += 100;
    l += 100;

    fill( k.lower_diag_rbegin(2), k.lower_diag_rend(2), 987.1203 );
    fill( k.col_rbegin(5), k.col_rend(5), 0.0039848 ); 
    fill( k.row_rbegin(5), k.row_rend(5), 9990.0039848 ); 

    cout << "k^(-1)=\n" << k.inverse() << endl;
    cout << "det(k)=\n" << k.det() << endl;
    cout << "l^(-1)=\n" << l.inverse() << endl;
    cout << "det(l)=\n" << l.det() << endl;
  
    dynamic_matrix<long double> o(k/l);

    cout << "k=\n" << k << endl;
    cout << "l=\n" << l << endl;
    cout << "o=k/l=\n " << k/l << endl;
    cout << "o*l = \n" << o*l << endl;

    dynamic_matrix<long double> m(17,17);
    fill( m.rbegin(), m.rend(), 990 );
    cout << "m=\n" << m << endl;
    cout << "m^(-1)=\n" << m.inverse() << endl;
    cout << "det(m)=\n" << m.det() << endl;

    dynamic_matrix<long double> n(2, 2);
    fill( n.rbegin(), n.rend(), 0 );
    fill( n.diag_rbegin(), n.diag_rend(), 100 );
    n[0][0]= -1;
    cout << "n = \n" << n << endl;
    cout << "det(n) = \n" << n.det() << endl;

    dynamic_matrix<long double> p( 101, 101 );
	feng::variate_generator<long double> vg;
	copy( vg.begin(), vg.begin()+101*101, p.begin());
    fill( p.rbegin(), p.rend(), vg() );
    fill( p.diag_rbegin(), p.diag_rend(), 190 );
    fill( p.upper_diag_rbegin(7), p.upper_diag_rend(7), 0.91 );
    fill( p.lower_diag_rbegin(4), p.lower_diag_rend(4), 1.9 );
    cout << "p=\n" << p << endl;
    cout << "p^(-1) = \n" << p.inverse() << endl;
    cout << "det(p) = \n" << p.det() << endl;
    dynamic_matrix<long double> p_inv = p.inverse();

    cout << setprecision(3);
    cout << "p*p^(-1) = \n" << p * p_inv << endl;
    
	dynamic_matrix<long double> q( 128, 128 );
	copy( vg.begin(), vg.begin()+128*128, q.begin());
    fill( q.diag_rbegin(), q.diag_rend(), 190 );
    fill( q.upper_diag_rbegin(7), q.upper_diag_rend(7), 0.1 );
    fill( q.lower_diag_rbegin(4), q.lower_diag_rend(4), 1.1 );
    cout << "q.\n" << q<< endl;
    cout << "q.(-1) = \n" << q.inverse() << endl;
    cout << "det(q. = \n" << q.det() << endl;


    dynamic_matrix<long double> q_inv = q.inverse();

    cout << setprecision(3);
    cout << "q*q^(-1) = \n" << q * q_inv << endl;

    return 0;
}





